Benjamin A. Levy1, Christopher M. Legault2, Timothy J. Miller2, Elizabeth N. Brooks2
1Ocean Associates Inc. in support of 2National Marine Fisheries Service, Northeast Fisheries Science Center, 166 Water Street, Woods Hole, Massachusetts 02543-1026
Corresponding author: Ben Levy (blevy6@gmail.com)
Competing interests: The authors declare there are no competing interests.
Many stock assessments in the United States use biomass estimates derived from stratified random bottom trawl data. To accurately track changes in true biomass of a species, catch observations must contain low noise and the proportion of the population that is sampled should be consistent over time. One increasingly documented phenomenon that violates the latter assumption is climate-induces movement leading to changes in relative biomass between sampled and unsampled portions of the stock area. Using the R package MixFishSim, we have developed data-driven spatial models for Yellowtail Flounder, Atlantic Cod, and Haddock in the western Atlantic Ocean to evaluate two methods of estimating population trend given temperature driven movement. Movement rates combine species-specific static habitat preferences with temperature tolerances. Habitat preferences were derived from niche models relating bottom trawl catches to environmental covariates. A repeating yearly temperature pattern produces repeating spatial biomass distributions in a given week, while a temperature gradient that increases on average over time results in spatial preferences that evolve throughout a given simulation. We simulated spatial time series data sets for each species for several temperature scenarios and population trends. Using stratified random sampling on model output we were able to compare biomass estimates derived from the design-based stratified mean to those using a spatio-temporal model-based approach that allows inclusion of environmental covariates. Our focus was on the ability of contemporary index standardization methods to track population trends under shifting spatial preferences.
survey, design-based, model-based, climate change, population simulation, VAST
The continental shelf off the coast of the Northeast United States spans from the Outer Banks of North Carolina to the Gulf of Maine. The region covers over 250,000 km\(^2\) of ocean, extending over 200 km from shore in the largest areas in New England to just 30 km off shore in the southern regions. This ecologically diverse region contains approximately 18,000 vertebrate marine species, making commercial fisheries an important part of local economies for centuries. In 2019, New England fisheries produced $22 billion in sales, which supported over 200,000 jobs (“Fisheries Economics of the United States: Data and Visualizations” n.d.). Maintaining a healthy ecosystem is therefore vital to the economic prosperity of the region.
The Northeast Fisheries Science Center (NEFSC) has conducted a bottom trawl survey since 1963 to support assessment and management of the fish and invertebrate populations in the region (Azarovitz 1981; Politis et al. 2014). The survey uses a stratified random design where bottom trawl sampling takes place in predefined strata along the eastern continental shelf, where depth and latitude define the strata. The survey has created a rich time series data set with many uses including species-specific habitat identification, age and life history information, analysis of how environmental conditions influence species biomass, and estimating yearly species biomass trends to help inform stock assessments, and ultimately fishing quota limits. The survey takes place twice each year- once in the spring and again in the fall. Most spatial analyses and projections of future distributions typically assume a constant survey catchability and availability over time. For this reason, NOAA’s survey design includes sampling each stratum in approximately the same 6 week time period in each season.
One component of a stock assessment is an estimation of yearly biomass that is often derived from the bottom trawl survey, when applicable. Stock assessment scientists choose from a number of approaches to obtain biomass estimates ranging from traditional design-based estimates to model-based estimates that vary in complexity. Design-based estimators rely on the design of the sampling scheme with the underlying assumption that the data being collected is representative of the population of interest. These methods typically account for spatial variation in samples by stratifying the sampling area. While it is possible to include environmental influences with design-based methods (Deville and Särndal 1992), it is not a common practice. Model-based biomass estimates use statistical models to measure the relationship between response variables (such as presence or biomass) and predictor variables (such as environmental factors). Model-based estimators help account for complex relationships between variables and can help overcome problems related to the sampling process. Common model-based approaches include Generalized Linear Models (GLM), Generalized Additive Model (GAM), and General Linear Mixed Models (GLMM).
Due to a combination of climate change and shifts in circulation, the Northeast United States continental shelf has experienced rapid warming in recent decades. At the same time a number of recent studies have found a shift in spatial distributions of many species (Nye et al. 2009; Henderson et al. 2017; Kleisner et al. 2017). Since stock assessment models rely on accurate descriptions of population dynamics and contemporary patterns of spatial abundance, there is concern that rapid undocumented changes in spatial distributions of species will bias future stock assessments. More specifically, as fish populations shift their distributions over time, catchability and/or availability in the survey could change, altering the relationship between the index and the true population (Arreguı́n-Sánchez 1996; Langan et al. 2021). A species shifting its range beyond the survey area or outside of the assumed stock boundary is an additional compounding factor to consider. Existing research has focused on temperature as the driver of such changes (Klein et al. 2017) and evidence suggests that failing to account for the impact of climate-induced change can lead to management challenges (Kerr et al. 2022). In these scenarios, management strategy evaluations have shown that misconceptions of stock status can lead to unintended overfishing, which can ultimately have detrimental ecologic and economic impacts (Mazur et al. 2023). We are therefore interested in analyzing the impact of climate driven movement and distribution changes on the accuracy of biomass estimates derived from NOAA’s ongoing bottom-trawl survey along the east coast.
To test the ability of the bottom trawl survey to track population trends under shifting environmental conditions, we construct spatial models for fish where movement depends on temperature preferences. We consider the impact of climate change by comparing simulations that use a repeating water temperature pattern to those where temperature increases on average over time. In both cases we conduct stratified random sampling (SRS) on model output and analyze the ability of the samples to track population trends. We compare yearly biomass estimates obtained from the stratified mean to estimates obtained from the Vector-Autoregressive Spatio-temporal (VAST) model. The stratified mean is a design-based approach that calculates the stratum specific mean catch per tow, which are then weighted by stratum area to obtain an overall mean. VAST is a spatial delta-generalized linear mixed model that estimates both biomass and probability of occurrence (presence/absence) (Thorson 2019). VAST also allows users to include covariate data to better inform the model. Covariates can be static (e.g. habitat preferences) or dynamics (e.g. temperature). We explore whether including environmental predictors leads to improved biomass estimates, which is particularly relevant as climate change progresses.
We construct spatial population models for Yellowtail Flounder (Limanda ferruginea), Atlantic Cod (Gadus morhua), and Haddock (Melanogrammus aeglefinus) on Georges Bank, hereafter referred to as Yellowtail, Cod, and Haddock. Movement of each species combines static species-specific habitat preferences with biologically-based temperature preferences. Model dynamics are driven by a time series of temperature gradients to create simulated spatial data sets for each population where the true biomass is known. Using temperature gradients that repeat each year creates data sets with repeating spatial patterns, whereas using a temperature gradient that increases on average throughout the simulation leads to spatial distributions that shift over time. We conduct SRS on our simulation data to mimic the bottom trawl survey and use the samples to compare the ability of contemporary index standardization methods to track population trends.
We use the R package MixFishSim (MFS) to model our populations (Dolder et al. 2020). MFS is a discrete spatiotemporal simulation tool that allows users to model multiple species under varying environmental conditions. The package uses a delay-difference population model with discrete processes for growth, death, and recruitment of the population. To address our research question we specify the following inputs for the MFS package.
Study Area
A shapefile for the 15 strata that comprise Georges Bank was discretized into a raster with 88 rows and 144 columns to use as our modeling environment. Each cell in our simulation domain represents approximately 8.7 km\(^2\). A fish stock is considered to be a subpopulation of a species that has similar intrinsic parameters. Each of the species being modeled has multiple biologically distinct stocks along the Atlantic coast resulting from local environmental conditions. As a result, each stock inhabits a different number of strata on Georges Bank. Haddock inhabit all 15 strata in the domain, Cod populate 13 strata, and Yellowtail can be found in 9 strata (Figure @ref(fig:strata-plot)).
(ref:crossref) (Walsh et al. 2015)
Strata inhabited by each species in our population models. All strata on the eastern US contintental shelf are shown in the first image with Georges Bank highlighted in orange (ref:crossref). Stratum numbers used by the NEFSC bottom trawl survey are shown in the middle row. Strata that are excluded from certain stratified random sampling are shown in the bottom row in yellow.
Population Dynamics and Recruitment
The time step for our models is one week. MFS uses a modified two-stage Deriso-Schnute delay difference equation that models the biomass in each cell in our study area (Dolder et al. 2020). Parameters in the population model account for growth of mature adults, mortality (natural and fishing), and the addition of new recruits. Recruitment is a function of the adult biomass that existed in the previous year and is added to the population incrementally throughout each species’ predefined spawning period. Parameter inputs were either obtained from the literature or chosen to produce desired model dynamics, and therefore, these results cannot be compared directly with stock assessment results. A full list of parameters used in our population model can be seen in Tables 1 and @ref(tab:paramsSCENARIOS).
Movement
The package combines species-specific temperature tolerances with habitat preferences to drive the probability of movement from cell \(I\) to cell \(J\) in week \(wk+1\) using the formulation
\[\begin{align} Pr(C_{wk+1}=J|C_{wk}=I) = \frac{e^{-\lambda \cdot d_{I,J}}\cdot(Hab^2_{J,s} \cdot Tol_{J,s,wk})}{\sum^C_{c=1}e^{-\lambda \cdot d} \cdot (Hab^2_{c,s} \cdot Tol_{c,s,wk})}, \label{moveP} \end{align}\]
where
\(e^{-\lambda \cdot d_{I,J}}\) accounts for distance between cells \(I\) and \(J\),
\(Hab^2_{J,s}\) is the static habitat value for species \(s\) in cell \(J\), and
\(Tol_{c,s,wk}\) is the value from normally distributed temperature tolerance for species \(s\) in cell \(c\) in week \(wk\).
The MFS package was designed to generate hypothetical temperature gradients and theoretical habitat preferences using Gaussian Random Fields. The following sections describe how we configured the habitat and temperature components to model our three species on Georges Bank.
Habitat Input
The habitat input for each species in MFS should represent spatial habitat preferences for each species that are time-invariant and are instead related to static geographic features. Species-specific habitat preferences were derived from niche models for each species using the lrren tool from the R package envi (Buller 2022). The lrren tool estimates an ecological niche using the relative risk function by relating presence/absence data to two covariate predictors. We used bottom trawl point data from 2009-2021 as our presence/absence input by using a value of 0 for any tow that failed to catch the given species and by representing a successful catch by the biomass of the given tow. We combined data from both the fall and spring surveys to obscure seasonal influences to allow the niche model to instead infer static habitat preferences independent of temperature. Depth and mean sediment size were used as our covariate predictors. Estimated depth for the region was obtained from FVCOM (Chen et al. 2006). The mean sediment size raster was interpolated in ArcMap using the natural neighbor interpolation method using point data collected by the United States Geologic Survey (USGS) (McMullen et al. n.d.). Since the values in \(Hab_{J,s}\) are required to be between 0 and 1, we rescaled the spatial estimates from lrren to fall between these bounds. See Figure @ref(fig:hab-plot1) for a visual representation of this process being applied to Cod. Figure @ref(fig:hab-plot) depicts habitat preferences \(Hab_{J,s}\) for each species.
Visual representation of niche model for Atlantic Cod. In the left most plot, tows with positive Cod catch are shown in red while tows without Cod are in blue. Positive tows are weighted by the biomass of the given catch. The middle two plots show our covariate predictors used in the niche model- depth and median sediment size. The plot on the right shows the final niche model that has been rescaled between 0 and 1, where blue values are closer to 0 and red values are closer to 1.
Static habitat preferences for each species in our population models. From left to right: Yellotwtail Flounder, Atlantic Cod, and Haddock.
Temperature Input
Each species is assumed to have normally distributed temperature preferences \(N(\mu,\sigma)\), with mean \(\mu\) and standard deviation \(\sigma\). Temperature values are in degrees Celcius. Values were chosen by combining information in the literature with temperatures recorded in the bottom trawl survey. We assume Yellowtail’s preferences are \(N(8.75,4.25)\), while Haddock and Cod have preferences \(N(9,4)\). Weekly estimated temperature data for the region for 2012 was obtained from FVCOM (Chen et al. 2006). We chose to repeat temperature estimates for a single year rather than use data for consecutive years to reduce the number of factors impacting model dynamics while still incorporating real data. The 2012 data was chosen because it displayed an average temperature pattern that consistently oscillated between maximum and minimum temperature values, allowing for a smooth repeating yearly temperature pattern for the constant temperature scenario. The 2012 temperature data was also transformed to create an oscillating pattern that increases 5 degrees Celsius on average over the duration of the simulation. We chose a 5 degree increase over a 20 year simulation to allow temperature change to have a meaningful impact on spatial patterns while remaining within reasonable computational limits in terms of the length of the simulation. Figure @ref(fig:temp-scenarios) depicts mean trends for the temperature scenarios used in our population models. dont forget to include animated gif in final submission